Ειδικά Μαθήματα Βοτανικής

logo

Όλα τα απαιτούμενα δεδομένα και αρχεία θα τα βρείτε στον ιστότοπο του μαθήματος στο e-class

1 Προετοιμασία

1.1 Εγκατάσταση βιβλιοθήκης

Εγκατάσταση του πρώτου και κύριου πακέτου (και ενός βοηθητικού)
(Αυτό γίνεται μια φορά σε κάθε Η/Υ - εάν έχει εγκατασταθεί δεν χρειάζεται να δοθεί ξανά η εντολή αυτή)

## CRAN repository:
install.packages("easypackages", repos = "http://cran.us.r-project.org")

## Github:
devtools::install_github("drsimonj/twidlr")

1.2 Φόρτωση βιβλιοθήκης

Φόρτωση του εν λόγω πακέτου

library(easypackages)

1.3 Εγκατάσταση βιβλιοθηκών

Με την εντολή packages('insert_package_name_here') γίνεται η εγκατάσταση των πακέτων που μας ενδιαφέρουν1

Δοκιμάστε το

1.4 Φόρτωση βιβλιοθηκών

Εντολή για φόρτωση των πακέτων

libraries("car", "psych", "tidyverse", "broom", "twidlr", "leaps", "relaimpo", 
    "gvlma", "MASS")

1.5 Working Directory

Με την ακόλουθη εντολή, βλέπουμε σε ποιόν φάκελο δουλεύουμε:

getwd()

1.6 Αλλαγή του working directory

Ενώ με την εντολή αυτή, αλλάζουμε τον φάκελο στον οποίο δουλεύουμε και αποθηκεύουμε δεδομένα
[π.χ., setwd("E:/") εάν θέλουμε να δουλεύουμε στον σκληρό δίσκο Ε]

setwd("E:/Academic/Lessons/<c5><e9><e4><e9><ea><dc> <cc><e1><e8><de><ec><e1><f4><e1> <c2><ef><f4><e1><ed><e9><ea><de><f2>/Labs/Labs")

2 Δεδομένα προς ανάλυση

2.1 Εισαγωγή αρχείου

Τώρα ας εισάγουμε το αρχείο το οποίο περιέχει τα δεδομένα τα οποία μας ενδιαφέρουν να αναλύσουμε

Aegean <- readxl::read_excel("G:/Academic/Lessons/EMB/Labs/Labs/Aegean Islands.xlsx")

Κάθε αρχείο πρέπει να έχει ένα όνομα στην R
Ποιό είναι το όνομα του αρχείου μας;

2.2 Δομή των δεδομένων μας

Ας δούμε τα δεδομένα τα οποία εισάγαμε στην R

Aegean

και τι είδους δεδομένα είναι αυτά (ποιοτικά; ποσοτικά;)

str(Aegean)
## Classes 'tbl_df', 'tbl' and 'data.frame':    59 obs. of  14 variables:
##  $ Island   : chr  "Alonissos" "Amorgos" "Anafi" "Andros" ...
##  $ Area     : num  65 121 38 380 19.7 ...
##  $ Elevation: num  475 823 584 1003 378 ...
##  $ MAT      : num  111 133 139 130 139 ...
##  $ MAP      : num  489 498 498 393 362 ...
##  $ Temp     : num  161 175 181 170 181 ...
##  $ Rain     : num  495 528 523 463 402 ...
##  $ Dis      : num  4.76 10.33 21.14 9.29 1.7 ...
##  $ Dm       : num  41 117.8 207.1 55.6 112.2 ...
##  $ Geology  : num  5 5 8 6 3 1 5 4 1 10 ...
##  $ Total    : num  515 728 504 915 311 ...
##  $ Regional : num  21 75 45 58 28 8 24 52 34 96 ...
##  $ Endemics : num  18 59 35 47 25 6 15 36 17 30 ...
##  $ SIE      : num  0 1 0 3 1 0 0 1 0 1 ...

Μπορούσαμε να πάρουμε την πληροφορία αυτή με άλλον τρόπο;

2.3 Λιγότερα δεκαδικά

Με την ακόλουθη εντολή, θα εμφανίζονται πλέον λιγότερα δεκαδικά:

options(digits = 2)

3 Προκαταρκτικές αναλύσεις

3.1 Έλεγχος κανονικότητας

3.1.1 Α

Δημιουργούμε ένα αρχείο το οποίο ελέγχει εάν όλες οι μεταβλητές μας (εκτός της πρώτης στήλης, η οποία είναι ποσοτική μεταβλητή - το όνομα των νησιών), έχουν κανονική κατανομή:

normality <- lapply(Aegean[, 2:14], shapiro.test)

3.1.2 Β

Ας δούμε τα αποτελέσματα του προηγούμενου βήματος

normality
## $Area
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.3, p-value = 3e-15
## 
## 
## $Elevation
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 2e-06
## 
## 
## $MAT
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 1e-07
## 
## 
## $MAP
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 2e-04
## 
## 
## $Temp
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 3e-07
## 
## 
## $Rain
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 5e-05
## 
## 
## $Dis
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.8, p-value = 5e-07
## 
## 
## $Dm
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 0.01
## 
## 
## $Geology
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 1e-04
## 
## 
## $Total
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.9, p-value = 3e-04
## 
## 
## $Regional
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.6, p-value = 2e-11
## 
## 
## $Endemics
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.5, p-value = 3e-13
## 
## 
## $SIE
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.2, p-value = 3e-16

Οι μεταβλητές με πιθανότητα μικρότερη από 0.05 ΔΕΝ εμφανίζουν κανονική κατανομή.

3.2 Έλεγχος συσχέτισης

Ας ελέγξουμε τη συσχέτιση μεταξύ των μεταβλητών μας
Εάν κάποιες έχουν συσχέτιση \(\geq 0.8-0.9\), πρέπει να τις αφαιρέσουμε από τη μεταγενέστερη ανάλυση (Γιατί;)

corr.test(Aegean[, 2:14], method = "spearman")
## Call:corr.test(x = Aegean[, 2:14], method = "spearman")
## Correlation matrix 
##            Area Elevation   MAT   MAP  Temp  Rain   Dis    Dm Geology
## Area       1.00      0.75 -0.39  0.10 -0.42  0.07  0.15 -0.41    0.79
## Elevation  0.75      1.00 -0.33  0.07 -0.38  0.04  0.04 -0.32    0.67
## MAT       -0.39     -0.33  1.00  0.00  0.96  0.05 -0.23  0.46   -0.32
## MAP        0.10      0.07  0.00  1.00  0.05  0.99  0.40 -0.26    0.03
## Temp      -0.42     -0.38  0.96  0.05  1.00  0.09 -0.19  0.33   -0.36
## Rain       0.07      0.04  0.05  0.99  0.09  1.00  0.38 -0.22    0.01
## Dis        0.15      0.04 -0.23  0.40 -0.19  0.38  1.00 -0.27    0.22
## Dm        -0.41     -0.32  0.46 -0.26  0.33 -0.22 -0.27  1.00   -0.23
## Geology    0.79      0.67 -0.32  0.03 -0.36  0.01  0.22 -0.23    1.00
## Total      0.95      0.77 -0.42  0.09 -0.43  0.05  0.16 -0.47    0.79
## Regional   0.63      0.65 -0.19  0.11 -0.24  0.10  0.10 -0.17    0.58
## Endemics   0.43      0.45 -0.03 -0.15 -0.12 -0.15 -0.07  0.18    0.45
## SIE        0.62      0.53 -0.31  0.11 -0.32  0.09  0.35 -0.32    0.50
##           Total Regional Endemics   SIE
## Area       0.95     0.63     0.43  0.62
## Elevation  0.77     0.65     0.45  0.53
## MAT       -0.42    -0.19    -0.03 -0.31
## MAP        0.09     0.11    -0.15  0.11
## Temp      -0.43    -0.24    -0.12 -0.32
## Rain       0.05     0.10    -0.15  0.09
## Dis        0.16     0.10    -0.07  0.35
## Dm        -0.47    -0.17     0.18 -0.32
## Geology    0.79     0.58     0.45  0.50
## Total      1.00     0.67     0.48  0.67
## Regional   0.67     1.00     0.77  0.53
## Endemics   0.48     0.77     1.00  0.49
## SIE        0.67     0.53     0.49  1.00
## Sample Size 
## [1] 59
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##           Area Elevation  MAT  MAP Temp Rain  Dis   Dm Geology Total
## Area      0.00      0.00 0.11 1.00 0.06 1.00 1.00 0.06    0.00  0.00
## Elevation 0.00      0.00 0.47 1.00 0.14 1.00 1.00 0.53    0.00  0.00
## MAT       0.00      0.01 0.00 1.00 0.00 1.00 1.00 0.02    0.54  0.06
## MAP       0.46      0.61 0.98 0.00 1.00 0.00 0.08 1.00    1.00  1.00
## Temp      0.00      0.00 0.00 0.68 0.00 1.00 1.00 0.47    0.25  0.04
## Rain      0.58      0.77 0.73 0.00 0.50 0.00 0.14 1.00    1.00  1.00
## Dis       0.27      0.75 0.08 0.00 0.15 0.00 0.00 1.00    1.00  1.00
## Dm        0.00      0.01 0.00 0.05 0.01 0.09 0.04 0.00    1.00  0.01
## Geology   0.00      0.00 0.01 0.84 0.01 0.95 0.10 0.08    0.00  0.00
## Total     0.00      0.00 0.00 0.50 0.00 0.69 0.22 0.00    0.00  0.00
## Regional  0.00      0.00 0.14 0.39 0.07 0.45 0.45 0.20    0.00  0.00
## Endemics  0.00      0.00 0.79 0.26 0.36 0.25 0.62 0.17    0.00  0.00
## SIE       0.00      0.00 0.02 0.39 0.01 0.50 0.01 0.01    0.00  0.00
##           Regional Endemics  SIE
## Area             0     0.04 0.00
## Elevation        0     0.02 0.00
## MAT              1     1.00 0.62
## MAP              1     1.00 1.00
## Temp             1     1.00 0.53
## Rain             1     1.00 1.00
## Dis              1     1.00 0.27
## Dm               1     1.00 0.53
## Geology          0     0.02 0.00
## Total            0     0.01 0.00
## Regional         0     0.00 0.00
## Endemics         0     0.00 0.01
## SIE              0     0.00 0.00
## 
##  To see confidence intervals of the correlations, print with the short=FALSE option

Διακρίνετε κάποιες μεταβλητές με τόσο μεγάλη συσχέτιση;

3.2.1 Αυτοματοποιημένος τρόπος

usdm::vifcor(Aegean %>% as.data.frame() %>% dplyr::select(Area:Geology), th = 0.8)
## 2 variables from the 9 input variables have collinearity problem: 
##  
## MAT MAP 
## 
## After excluding the collinear variables, the linear correlation coefficients ranges between: 
## min correlation ( Rain ~ Elevation ):  -0.021 
## max correlation ( Geology ~ Elevation ):  0.75 
## 
## ---------- VIFs of the remained variables -------- 
##   Variables VIF
## 1      Area 2.6
## 2 Elevation 3.9
## 3      Temp 1.9
## 4      Rain 1.3
## 5       Dis 1.7
## 6        Dm 1.2
## 7   Geology 2.4

3.3 Οπτικοποίηση αποτελεσμάτων

Μπορούμε να οπτικοποιήσουμε τα αποτελέσματα μας και να διαπιστώσουμε εάν κάποια μεταβλητή δεν εμφανίζει κανονική κατανομή

scatterplotMatrix(Aegean[, 2:14], spread = F, smoother.args = list(lty = 2), 
    main = "Scatter Plot Matrix")

3.4 Κανονικοποίηση των δεδομένων μας

Καθώς τα δεδομένα μας δεν έχουν κανονική κατανομή, θα πρέπει να τα κανονικοποιήσουμε. Αυτό σημαίνει ότι θα πρέπει να τα λογαριθμήσουμε. Συνεπώς, πρέπει να ελέγξουμε εάν κάποιες μεταβλητές περιέχουν το 0. Η ακόλουθη εντολή μας βοηθάει να το διαπιστώσουμε αυτό:

range(Aegean$SIE)
## [1]   0 172

Κάντε το και για τις υπόλοιπες μεταβλητές (Δηλαδή: range(Aegean$insert_IV_name_here))

3.5 Λογαρίθμηση

Εφ’όσον έχουμε ελέγξει όλες τις μεταβλητές μας, μπορούμε να τις λογαριθμήσουμε:

Aegean <- Aegean %>% mutate(Area = log10(Area), Total = log10(Total), Endemics = log10(Endemics), 
    Regional = log10(Regional), SIE = log10(SIE + 1))

Και εάν θέλουμε να λογαριθμήσουμε όλες τις μεταβλητές που μας ενδιαφέρουν, ένας γρήγορος τρόπος είναι ο ακόλουθος:

Aegean_log <- lapply(Aegean[, 2:13], as.numeric) %>% as_tibble() %>% log10() %>% 
    dplyr::mutate(Island = Aegean$Island, SIE = log10(Aegean$SIE + 1)) %>% dplyr::select(Island, 
    everything(), SIE)

3.6 Δομή των δεδομένων μας (ξανά)

Ας δούμε τα δεδομένα μας ξανά:

Aegean_log

4 Κυρίως ανάλυση

4.1 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο enter

Τώρα είμαστε πλέον σε θέση να δημιουργήσουμε 4 μοντέλα πολλαπλής γραμμικής παλινδρόμησης, ένα για κάθε κατηγορία φυτικών taxa [Συνολικός αριθμός ειδών, Ενδημικά είδη, Περιφερειακά ενδημικά είδη και Τοπικά νησιωτικά είδη (Total, Endemics, Regional & SIE, αντίστοιχα)] σε σχέση με όλες τις ανεξάρτητες μεταβλητές μας2. Θα χρησιμοποιήσουμε αρχικά την μέθοδο enter:

enter1 <- lm(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, data = Aegean_log)

enter2 <- lm(Endemics ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log)

enter3 <- lm(SIE ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, data = Aegean_log)

enter4 <- lm(Regional ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log)

ΠΡΟΣΟΧΗ
Τα ανωτέρω μοντέλα δεν δείχνουν ποιές είναι οι μεταβλητές οι οποίες επηρεάζουν τα πρότυπα ποικιλότητας. Δείχνουν απλά τι ποσοστό εξηγείται εάν συμπεριλαμβάνονται όλες οι ανεξάρτητες μεταβλητές στο μοντέλο μας

Ας δούμε τώρα την περίληψη για κάθε ένα εξ αυτών των μοντέλων:

summary(enter1)
## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30188 -0.03583 -0.00143  0.04867  0.16530 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.73748    1.01870    1.71    0.094 .  
## Area         0.23040    0.03857    5.97  2.3e-07 ***
## Elevation    0.01921    0.07729    0.25    0.805    
## Rain         0.12354    0.11986    1.03    0.308    
## Temp         0.06476    0.41385    0.16    0.876    
## Dm          -0.03509    0.02071   -1.69    0.096 .  
## Dis          0.00981    0.02605    0.38    0.708    
## Geology      0.15085    0.06651    2.27    0.028 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.088 on 51 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.834 
## F-statistic: 42.7 on 7 and 51 DF,  p-value: <2e-16
summary(enter2)
## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.501 -0.130  0.030  0.135  0.554 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -1.7112     3.0680   -0.56     0.58
## Area          0.1802     0.1161    1.55     0.13
## Elevation     0.3017     0.2328    1.30     0.20
## Rain         -0.3774     0.3610   -1.05     0.30
## Temp          1.2387     1.2464    0.99     0.33
## Dm            0.0821     0.0624    1.32     0.19
## Dis          -0.0205     0.0785   -0.26     0.79
## Geology       0.1331     0.2003    0.66     0.51
## 
## Residual standard error: 0.26 on 51 degrees of freedom
## Multiple R-squared:  0.383,  Adjusted R-squared:  0.298 
## F-statistic: 4.52 on 7 and 51 DF,  p-value: 0.000558
summary(enter3)
## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5765 -0.1688 -0.0096  0.1368  0.9108 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -0.4846     3.4535   -0.14   0.8890   
## Area          0.3987     0.1307    3.05   0.0036 **
## Elevation     0.4987     0.2620    1.90   0.0627 . 
## Rain         -0.1677     0.4063   -0.41   0.6816   
## Temp         -0.3976     1.4030   -0.28   0.7780   
## Dm           -0.0172     0.0702   -0.25   0.8070   
## Dis           0.1056     0.0883    1.20   0.2375   
## Geology      -0.1697     0.2255   -0.75   0.4553   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3 on 51 degrees of freedom
## Multiple R-squared:  0.602,  Adjusted R-squared:  0.547 
## F-statistic:   11 on 7 and 51 DF,  p-value: 2.18e-08
summary(enter4)
## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.5901 -0.0956  0.0419  0.1268  0.4315 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -3.27429    2.48522   -1.32    0.194  
## Area         0.23149    0.09408    2.46    0.017 *
## Elevation    0.31825    0.18856    1.69    0.098 .
## Rain         0.22590    0.29241    0.77    0.443  
## Temp         1.25730    1.00964    1.25    0.219  
## Dm           0.02453    0.05053    0.49    0.629  
## Dis          0.00645    0.06356    0.10    0.920  
## Geology      0.10430    0.16226    0.64    0.523  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.21 on 51 degrees of freedom
## Multiple R-squared:  0.557,  Adjusted R-squared:  0.496 
## F-statistic: 9.15 on 7 and 51 DF,  p-value: 2.83e-07

4.2 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο της βηματικής πολλαπλής γραμμικής παλινδρόμησης

Με την μέθοδο αυτή θα μπορέσουμε να δούμε ποιές μεταβλητές εμπερικλείονται στο βέλτιστο μοντέλο για τον συνολικό αριθμό των ειδών

stepm1 <- stepAIC(lm(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log), scale = 0, direction = "both", trace = 1, keep = NULL, 
    steps = 1000, k = 2)

summary(stepm1)

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιές μεταβλητές εξηγούν και σε τι ποσοστό τα πρότυπα της ποικιλότητας για αυτές τις κατηγορίες;

ΠΡΟΣΟΧΗ: Μας ενδιαφέρει το adjusted R squared

4.3 Ποιό μοντέλο είναι καλύτερο;

Με την βοήθεια του Akaike Information Criterion (εν συντομία: AIC), μπορούμε να διαπιστώσουμε ποιό μοντέλο είναι “καλύτερο”. Όσο μικρότερη η τιμή του AIC, τόσο καλύτερο το μοντέλο.

AIC(enter1, stepm1)

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιό από τα δύο μοντέλα για κάθε κατηγορία ειδών είναι καλύτερο;

4.4 Δημιουργία μιας λειτουργίας για τον υπολογισμό της σχετικής σημαντικότητας των μεταβλητών

relweights <- function(fit, ...) {
    R <- cor(fit$model)
    nvar <- ncol(R)
    rxx <- R[2:nvar, 2:nvar]
    rxy <- R[2:nvar, 1]
    svd <- eigen(rxx)
    evec <- svd$vectors
    ev <- svd$values
    delta <- diag(sqrt(ev))
    lambda <- evec %*% delta %*% t(evec)
    lambdasq <- lambda^2
    beta <- solve(lambda) %*% rxy
    rsquare <- colSums(beta^2)
    rawwgt <- lambdasq %*% beta^2
    import <- (rawwgt/rsquare) * 100
    import <- as.data.frame(import)
    row.names(import) <- names(fit$model[2:nvar])
    names(import) <- "Weights"
    import <- import[order(import), 1, drop = FALSE]
    dotchart(import$Weights, labels = row.names(import), xlab = "% of R-Square", 
        pch = 19, main = "Relative Importance of Predictor Variables", sub = paste("Total R-Square=", 
            round(rsquare, digits = 3)), ...)
    return(import)
}

4.5 Υπολογισμός της σχετικής σημαντικότητας των μεταβλητών

Μετά την περισπωμένη (~) βάζουμε όποιες μεταβλητές έχουν βγει σημαντικές από την βηματική πολλαπλή γραμμική παλινδρόμηση:

rel.1 <- lm(Total ~ Area + Dm + Geology, data = Aegean_log)

relweights(rel.1, col = "blue")
Weights
Dm 11
Geology 35
Area 54

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιά μεταβλητή είναι σημαντικότερη για την εκάστοτε κατηγορία ειδών;

4.6 Πολλαπλή γραμμική παλινδρόμηση με την μέθοδο των ολικών υποσυνόλων

Μέχρι στιγμής, έχουμε δει δύο διαφορετικές μεθόδους, όσον αφορά την πολλαπλή γραμμική παλινδρόμηση. Όμως, δεν μπορούμε να είμαστε 100% σίγουροι ότι έχουμε φθάσει στο ορθό συμπέρασμα, καθότι ένα από τα μειονεκτήματα της βηματικής πολλαπλής παλινδρόμησης, είναι ότι μπορεί να “κολλήσει” σε ένα υποβέλτιστο πλατώ. Τον σκόπελο αυτό, μπορούμε να τον αποφύγουμε χρησιμοποιώντας την μέθοδο των ολικών υποσυνόλων:

bes.tot <- regsubsets(Total ~ Area + Elevation + Rain + Temp + Dm + Dis + Geology, 
    data = Aegean_log, nbest = 4)

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

4.7 Απεικόνιση των αποτελεσμάτων του προηγούμενου βηματος

Το γράφημα αυτό θα μας δείξει ποιός είναι ο πραγματικά βέλτιστος συνδυασμός των ανεξάρτητων μεταβλητών:

plot(bes.tot, scale = "adjr2")

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

tot <- lm(Total ~ Area + Rain + Dm + Geology, data = Aegean_log)

summary(tot)
## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30141 -0.03724  0.00481  0.04501  0.16857 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9077     0.3196    5.97  1.9e-07 ***
## Area          0.2345     0.0301    7.79  2.1e-10 ***
## Rain          0.1341     0.1126    1.19    0.239    
## Dm           -0.0364     0.0197   -1.85    0.069 .  
## Geology       0.1543     0.0636    2.43    0.019 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.085 on 54 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.843 
## F-statistic: 78.9 on 4 and 54 DF,  p-value: <2e-16
relweights(tot, col = "blue")
Weights
Rain 0.62
Dm 10.45
Geology 35.16
Area 53.77

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional). Ποιές είναι οι σημαντικότερες μεταβλητές;

4.8 Έλεγχος των αποτελεσμάτων της παλινδρόμησης

Τώρα θα ελέγξουμε κατ’αρχάς εάν οι μεταβλητές του βέλτιστου μοντέλου μας παραβιάζουν μια από τις αρχές της πολλαπλής γραμμικής παλινδρόμησης:

sqrt(vif(tot)) > 2
##    Area    Rain      Dm Geology 
##   FALSE   FALSE   FALSE   FALSE

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

4.9 Έλεγχος πιθανοτήτων

Στο σημείο αυτό θα χρειαστεί να ελέγξουμε εάν η πιθανότητα κάθε μεταβλητής η οποία εμπερικλείεται στο τελικό μας μοντέλο είναι πραγματικά στατιστικώς σημαντική:

##--------------------------------------------------------------------------
## Let's have a look on the p-values for each of the IVs that have entered
## our final model
##--------------------------------------------------------------------------
summary(tot)
## 
## Call:
## stats::lm(formula = formula, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30141 -0.03724  0.00481  0.04501  0.16857 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9077     0.3196    5.97  1.9e-07 ***
## Area          0.2345     0.0301    7.79  2.1e-10 ***
## Rain          0.1341     0.1126    1.19    0.239    
## Dm           -0.0364     0.0197   -1.85    0.069 .  
## Geology       0.1543     0.0636    2.43    0.019 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.085 on 54 degrees of freedom
## Multiple R-squared:  0.854,  Adjusted R-squared:  0.843 
## F-statistic: 78.9 on 4 and 54 DF,  p-value: <2e-16
##--------------------------------------------------------------------------

##--------------------------------------------------------------------------
## Inside the parenthesis we include the above p-values
##--------------------------------------------------------------------------
p <- c(0, 0.239, 0.069, 0.019)
##--------------------------------------------------------------------------


##--------------------------------------------------------------------------
## Level of significance
##--------------------------------------------------------------------------
alp <- 0.05
##--------------------------------------------------------------------------


##--------------------------------------------------------------------------
## If the model's p-values are equal or smaller than the p-values we obtain
## from the following command, then they are ok
##--------------------------------------------------------------------------
p.adjust(p, method = "bonferroni", n = length(p))
## [1] 0.000 0.956 0.276 0.076
##--------------------------------------------------------------------------

Κάντε το ίδιο και για τις υπόλοιπες τρεις κατηγορίες ειδών (Endemics, SIE, Regional).

5 Εργασία για το σπίτι

Έχετε διορία επτά (7) ημερών να στείλετε την εργασία σας σε μορφή PDF3 σε αυτό το e-mail και να απαντήσετε στα ακόλουθα ερωτήματα, αφού έχετε περιγράψει το σετ δεδομενων το οποίο έχετε χρησιμοποιήσει:
1. Πόσα νησιά έχει;
2. Πού βρίσκονται;
3. Ποιο είναι το μεγαλύτερο και ποιο το μικρότερο;
4. Ποιο είναι το ψηλότερο και ποιο το χαμηλότερο;
5. Ποιο είναι πιο κοντά στην ηπειρωτική ξηρά και ποιο είναι πιο μακρυά;
6. Ποιο έχει τα περισσότερα και ποιο τα λιγότερα είδη;

Με βάση τα αποτελέσματα σας, απαντήστε στα ακόλουθα ερωτήματα:
7. Ποιές μεταβλητές δεν εμφανίζουν κανονική κατανομή4;
8. Ποιές ενέργειες θα κάνετε προκειμένου να κανονικοποιήσετε την κατανομή των μεταβλητών αυτών;
9. Ποιές μεταβλητές εμφανίζουν μεγάλο (\(>0.8\)) συντελεστή συσχέτισης;
10. Εάν υπάρχουν κάποια ζεύγη ανεξάρτητων μεταβλητών με σ.σ. \(>0.8\), θα τις χρησιμοποιήσετε όλες; 11. Για κάθε έναν από τους τύπους μοντέλων [Enter, βηματική πολλαπλή γραμμική παλινδρόμηση και ολικά υποσύνολα], καθώς και για όλες τις κατηγορίες ειδών (Total, Endemics, SIE, Regional), να αναφέρετε το adj. r2 και τις μεταβλητές που περιέχονται στην περίληψη των μοντέλων5.
12.Συμφωνούν τα αποτελέσματα μεταξύ της βηματικής πολλαπλής γραμμικής παλινδρόμησης και της μεθόδου των ολικών υποσυνόλων; Εάν όχι, που διαφέρουν;
13. Ποιο είναι το ποσοστό συμμετοχής των μεταβλητών (για κάθε κατηγορία ειδών) στα μοντέλα της βηματικής πολλαπλής γραμμικής παλινδρόμησης και της μεθόδου των ολικών υποσυνόλων;
14. Οι τιμές των πιθανοτήτων των μεταβλητών που βρίσκονται εντός του βέλτιστου μοντέλου, είναι ίσες ή μικρότερες από αυτές του Bonferroni correction;
15. Να φτιάξετε την εξίσωση για κάθε κατηγορία ειδών, όπως αυτή προκύπτει από το βέλτιστο μοντέλο.

Τα σετ δεδομένων για την εργαστηριακή άσκηση αυτή μπορείτε να τα βρείτε εδώ.


  1. Βλέπετε ποιά πακέτα μας ενδιαφέρουν από το επόμενο βήμα

  2. Ποιές ονομάζουμε ανεξάρτητες μεταβλητές; Μπορούμε να τις χρησιμοποιήσουμε όλες; Εάν όχι, γιατί;

  3. Διαφορετικά δεν θα γίνει δεκτή και δεν θα βαθμολογηθείτε

  4. Να συμπεριλάβετε και το αντίστοιχο γράφημα

  5. Τα αποτελέσματα να φαίνονται σε πίνακα τον οποίο θα έχετε διαμορφώσει κατάλληλα

Κώστας Κουγιουμουτζής

Spring semester 2019-2020

 

Delivered to you by Dr. Kostas Kougioumoutzis

kkougiou@aua.gr